	
\documentclass{article} % 文档类别: article, report, book, letter, 等等
\usepackage{ctex}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{listings}
\usetikzlibrary{intersections}
\usetikzlibrary{calc}
\title{刘徽割圆术将π的求解转为开方公式}
\author{朱正路}
\date{\today} % 使用当前日期，也可以指定特定日期
% 设置 listings 的样式，去掉语法高亮
\lstset{
	language=Python,
	basicstyle=\ttfamily, % 使用等宽字体
	backgroundcolor=\color{white}, % 背景色设为白色，去除阴影效果
	keywordstyle=\relax, % 关闭关键字高亮
	commentstyle=\relax, % 关闭注释高亮
	stringstyle=\relax, % 关闭字符串高亮
	showstringspaces=false, % 不在字符串中显示空格
	numbers=none, % 不显示行号
	breaklines=true, % 允许折行
	breakatwhitespace=true, % 在空白处折行
	tabsize=4, % 设置 tab 宽度
}
\begin{document}

	\maketitle % 创建标题页
	
	\section{引言}
	
	圆周率的计算，对于天文历法、气象都有重要的指导作用。刘徽割圆术利用圆的内接正多边形面积和来逐渐接近π值，但是王能超的《刘徽数学割圆术》，对于刘徽割圆术并没有一个递归的公式给出来，只有一个针对n多边形的实例计算。本文用公式表现抽象的数学，代替传统的实例。给出通用的基于平方根的递归式刘徽公式。
	\section{刘徽公式的计算方法}
	
	利用圆内接正n边形，定义该正n边形的每个三角形面积为\( \mathcal{S}_{\Delta AOC} \) 。
	为使用归纳法求解更精确的π，将内接正n边形变密1倍，定义变稠密后的正n边形的每个三角形面积为\( \mathcal{S}_{\Delta AOD} \)。
	
	利用 \( \mathcal{S}_{\Delta AOC} \) ，求解\( \mathcal{S}_{\Delta AOD} \)。这样已知正n多边形每个三角形的面积，就可以按归纳法求解正2n边形的每个三角形的面积。这里三角形是指以原点O构建的三角形，即图中的AOD和AOC三角形。

	假设圆为单位圆，即半径为1.
	图中角度ABO为直角。

	图中圆的面积为，
	\begin{equation}
	S = \pi \cdot r^2 = \pi     
\end{equation}

	图中正n边形的的面积为，
	\begin{equation}
S = n \cdot \mathcal{S}_{\Delta AOC} 
\end{equation}	

加密1倍后，图中正n边形的的面积为，
\begin{equation}
S = n \cdot \mathcal{S}_{\Delta AOD} 
\end{equation}
	
\begin{tikzpicture}[scale=2]
	% 重新定义圆的半径，放大1.2倍
	\def\radius{1.2}
	% 绘制放大后的圆并标记原点O
	\draw (0,0) coordinate[label=below left:$O$] (O) circle (\radius);
	
	% 定义点A、D、C，基于新的半径
	\coordinate (A) at (135:\radius); % 圆上点A，左上方向
	\coordinate (D) at (90:\radius);  % 圆上点D，正上方
	\coordinate (C) at (45:\radius);  % 圆上点C，右上方向
	
	% 绘制点A、D、C并标记
	\fill (A) circle (1.5pt) node[above left] {$A$};
	\fill (D) circle (1.5pt) node[above] {$D$};
	\fill (C) circle (1.5pt) node[above right] {$C$};
	
	% 绘制线段AO、OC、AD、DC、AC
	\draw (O) -- (A);
	\draw (O) -- (D);
	\draw (O) -- (C);
	\draw (A) -- (D);
	\draw (D) -- (C);
	\draw (A) -- (C);
	
	% 绘制直径AC，并找到其与OD的交点B
	\coordinate[label=below right:$B$] (B) at (intersection of O--D and A--C);
	\fill (B) circle (1.5pt);
	
	% 绘制线段OD和DB
	\draw (O) -- (D);
	\draw (D) -- (B);
	
	% 标记角ABO为直角，使用tkz-euclide包来更精确地控制角度标记
    % 手动添加直角标记
    \coordinate (ABMid) at ($(A)!0.9!(B)$); % AB的中点用于定位直角标记
    \coordinate (E) at ($(B) - (0,1mm)$);
	\draw ($(ABMid) - (0,1mm)$) -- ($(ABMid)$);
	\draw (E) -- ++(-0.8mm,0);
\end{tikzpicture}	
	\subsection{利用 \( \mathcal{S}_{\Delta AOC} \) ，求解\( \mathcal{S}_{\Delta AOD} \)。}
	求解稠密正多边形的每个三角形的面积
		\begin{equation}
		 \mathcal{S}_{\Delta AOD}  = \frac {1} {2} \cdot \overline{AB} \cdot \overline{OD} 
		  = \frac {1}{2} \cdot \overline{AB}
	\end{equation}
	这个容易看出来，由于∠ABO是直角，所以三角形AOD的面积如此。	
	从而建立了面积和AB线段长度的关系。
	
	求解正多边形的每个三角形的面积
		\begin{equation}
	\mathcal{S}_{\Delta AOC}  = \overline{AB} \cdot \overline{BO} 
\end{equation}

从而得到，
		\begin{equation}
	\overline{BO} = \frac {\mathcal{S}_{\Delta AOC}} {\overline{AB}} 
\end{equation}
	
	图中可以看到，
	
			\begin{equation}
		\overline{AB}^2 + \overline{BO}^2 =  \overline{AO}^2 = 1
	\end{equation}
	
	将公式（6）代入公式（7）
	
	从而得到，
	
	\begin{align}
		&\overline{AB}^2 + \frac {\mathcal{S}_{\Delta AOC} ^2} {\overline{AB}^2}  = 1 \\
		&\overline{AB}^4 - \overline{AB}^2 + \mathcal{S}_{\Delta AOC} ^2 = 0
	\end{align}
	
	利用韦达公式，可求解 AB线段长度。即
			\begin{equation}
	\overline{AB}^2 = \frac {1 \pm \sqrt {1 - 4 \cdot \mathcal{S}_{\Delta AOC} ^2 }} {2}
\end{equation}	

从而，

			\begin{equation}
	\overline{AB} = \sqrt {\frac {1 \pm \sqrt {1 - 4 \cdot \mathcal{S}_{\Delta AOC} ^2 }} {2} }
\end{equation}	

这里 AB线段长度有两个解，一大一小。从上面圆的标注图来看，分别对应 \(\overline {AB}\) 和 \(\overline {BO}\) 这两个解。
随着正n边形越来越稠密，\(\overline {AB}\)会越来越小于 \(\overline {BO}\)，故而\(\overline {AB}\)的长度取小值作为解。

从而，

			\begin{equation}
	\overline{AB} = \sqrt {\frac {1 - \sqrt {1 - 4 \cdot \mathcal{S}_{\Delta AOC} ^2 }} {2} }
\end{equation}	


将公式（4）代入，从而获得，

			\begin{equation}
	\mathcal{S}_{\Delta AOD} = \frac {1}{2} \cdot \sqrt {\frac {1 - \sqrt {1 - 4 \cdot \mathcal{S}_{\Delta AOC} ^2 }} {2} }
\end{equation}	

	这个公式表明，更小的正多边形原点三角形的面积可以由大的正多边形的原点三角形的面积，通过多次求平方根的方式获得。
	

为方便每次迭代，将下面公式称为刘徽公式。
			\begin{equation}
	\mathcal{S}_{2n} = \frac {1}{2} \cdot \sqrt {\frac {1 - \sqrt {1 - 4 \cdot \mathcal{S}_{n} ^2 }} {2} }
\end{equation}	

其中 \(\mathcal{S}_{n}\)为圆的内接正n边形的原点三角形的面积。
其中 \(\mathcal{S}_{2n}\)为圆的内接正2n边形的原点三角形的面积。

	联立公式（1）（3）（13）
可得到\(\pi\)的结果。
此时 \(\pi\)计算公式为，
	\begin{equation}
	\pi = n \cdot 	\mathcal{S}_{n}
\end{equation}	

为方便，令
			\begin{equation}
	\mathcal{B}_{n} = 1 -  4 \cdot \mathcal{S}_{n} ^2 
\end{equation}

则刘徽公式变为，

			\begin{equation}
	\mathcal{B}_{2n} = \frac {1} {2} + \frac {1} {2}  \cdot \sqrt {\mathcal{B}_{n} } 
\end{equation}

从而，

			\begin{equation}
	\mathcal{B}_{2n} = \frac {1} {2} + \frac {1} {2}  \cdot \sqrt {\frac {1} {2} + \frac {1} {2}  \cdot \sqrt {\frac {1} {2} + \frac {1} {2}  \cdot \sqrt {...}  }  } 	
\end{equation}

这里，
 			\begin{equation}
 	\mathcal{B}_{12} =  \frac {3} {4}
 \end{equation}
 
 定义如下方差，
 
  			\begin{equation}
 	\mathcal{\delta}_{4n}^2 =  \mathcal{\pi}_{4n}^2 - \mathcal{\pi}_{2n}^2
 \end{equation}
 
 将公式15,16,17代入上面得到，
 
   			\begin{equation}
 	\mathcal{\delta}_{4n} = 2n \cdot (1 - \mathcal{B}_{4n})
 \end{equation}
 
 这个方差可用于评估每次从正n边形稠密为正2n边形的时候计算π值的误差缩小的程度。
 
 此时，
 
   			\begin{equation}
 	\mathcal{\pi}_{n}^2 =  n^2 \cdot \frac {1- \mathcal{B}_{n}}{4}
 \end{equation}
 
 或者
    			\begin{equation}
 	\mathcal{\pi}_{4n}^2 =  4 \cdot n^2 \cdot (1- \mathcal{B}_{4n}) \end{equation}
 
 所以，刘徽公式比较适合计算π的平方。
 
		\subsection{计算实例}
计算当
n = 6

\begin{equation}
	\mathcal{S}_{6} = \frac {\sqrt{3}}{4} 
	\end{equation}	
这个原点小三角的面积可自行计算。
此时根据公式（1）（3）
	计算当
n = 6

\begin{equation}
	\pi = 6 \cdot 	\mathcal{S}_{6} = 6 \cdot \frac {\sqrt{3}}{4} = 2.5
	\end{equation}	
计算当n=12
			\begin{equation}
	\mathcal{S}_{12} = \frac {1}{2} \cdot \sqrt {\frac {1 - \sqrt {1 - 4 \cdot \mathcal{S}_{6} ^2 }} {2} } = \frac {1}{4}
\end{equation}		
	
	此时，
	
	\begin{equation}
		\pi = 12 \cdot 	\mathcal{S}_{12} = 12 \cdot \frac {1}{4} = 3
	\end{equation}	
计算当n=24
			\begin{equation}
	\mathcal{S}_{24} = \frac {1}{2} \cdot \sqrt {\frac {1 - \sqrt {1 - 4 \cdot \mathcal{S}_{12} ^2 }} {2} } = \frac {\sqrt{2-\sqrt{3}}}{4}
\end{equation}	
		此时，
	
	\begin{equation}
		\pi = 24 \cdot 	\mathcal{S}_{24} = 24 \cdot \frac {\sqrt{2-\sqrt{3}}}{4} = 3.1
	\end{equation}
	由此到n = 48, 96, 192, 384, 768, 1536, 3072...
		\subsection{使用计算机计算刘徽公式}
可使用python语言计算上面的程序，并评估精度,对比拉马努金公式或其他公式。

此时,

计算初始值使用公式19.

每次迭代使用公式17.

获得结果使用公式22.

评估精度使用公式21 （注意这个精度还不是和真正的π值的精度对比，只是和下一个正n边形的步长差)

代码如下
\begin{verbatim}

import math
n = 12  # 初始值
times = 10  # 迭代次数，可以根据需要调整
b0 = 3/4
b1 = 1/2 + 1/2*math.sqrt(b0)

for _ in range(times):
n *= 2  # 每次迭代n的值翻倍
print(n)
b0 = b1
b1 = 1/2 + 1/2 *math.sqrt(b0)

delta = n*(1-b1)
pai = n*math.sqrt(1-b1)
print("delta (difference between the next pai^2, not delta to the real pai^2): ", delta)
print(" pai: ", pai)

结果如下，
n :  12288
delta (difference between the next pai^2, not delta to the real pai^2):  0.0008031904444578686
pai:  3.1415926186407894

结果如下，
n :  98304
delta (difference between the next pai^2, not delta to the real pai^2):  0.00010039880726253614
pai:  3.1415926453212157
\end{verbatim}


当n更大时候，受限于python的处理双精度数据的方式上面代码处理会有问题，需要修改代码了。
\end{document}

